Classical simulation of noninteracting-fermion quantum circuits 

Barbara M. Terhal and David P. DiVincenzo 
IBM T.J. Watson Research Center, P.O. Box 218, Yorktown Heights, NY 10598 USA 



Abstract 

We show that a class of quantum computations that was recently shown to be 
efficiently simulatable on a classical computer by Valiant ||] corresponds to a 
physical model of noninteracting fermions in one dimension. We give an alter- 
native proof of his result using the language of fermions and extend the result 
to noninteracting fermions with arbitrary pairwise interactions, where gates 
can be conditioned on outcomes of complete von Neumann measurements in 
the computational basis on other fermionic modes in the circuit. This last 
result is in remarkable contrast with the case of noninteracting bosons where 
universal quantum computation can be achieved by allowing gates to be con- 
ditioned on classical bits 
PACS Numbers: 03.67.Lx, 5.30.-d 
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I. INTRODUCTION 



To understand the power of a quantum computer, it is worthwhile to explore under what 
restrictions that power is weakened so as to make the computation efficiently simulatable 
with the use of a classical device. A nontrivial example is the Knill-Gottesman class of 
quantum computations which can be shown to be efficiently simulatable on a classical 
computer. In the quantum circuit we only allow 1-qubit Hadamard transformations, 1- 
qubit 7r/2 phase shifts, 1-qubit Pauli-rotations and 2-qubit CNOT gates and furthermore 
the final measurements are projections in the two eigenspaces of any sequence of Pauli-matrix 
observables. 

This question about restricted classes of quantum computation is also related to the 
question of universality of a quantum computation. What set of gates, or in more physical 
terms, what physical system can be used to implement universal quantum computation? 
Surprises have been found in this direction, for example it was shown that the 2-qubit 
exchange interaction is sufficient for universal computation and, for example, that universal 
computation can be achieved with a network of phase shifters, beamsplitters and photon 
counters, i.e. noninteracting bosons, where logical gates can be conditioned on previous 
measurement outcomes 0. 

In Ref. a new class of quantum computations is introduced that are shown to be effi- 
ciently simulatable on a classical device. The class includes a special set of unitary 2-qubit 
gates on nearest-neighbor qubits. In this paper we will analyze this class of gates and show 
that it maps onto a system of noninteracting fermions (i.e., associated with Hamiltonian 
interactions that are quadratic in fermion creation and annihilation operators) in one di- 
mension. The equivalence will enable us to give a straightforward derivation of the classical 
simulation, as well as extend the class of quantum computations to include (1) noninteracting 
fermions without nearest neighbor restrictions and (2) gates that are applied conditionally on 
measurement outcomes. In particular the second result, when compared with the universal 
computation by linear optics in Ref. , shows a fundamental difference between bosons and 
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fermions. One possible cause for the difference is that the bosonic modes, unhke fermionic 
modes, can contain more than one particle, a feature that is employed in the universality 
construction in Ref. 0. 

The class of gates that Valiant shows to be classically simulatable is larger than the 
class that we start from in Theorem ^ (see below); his class includes non- unitary gates and 
also some special set of 2-qubit (possibly non-unitary) gates on the first 2 qubits. E. Knill 
has now shown that this entire class is (indeed) weaker than full quantum computation. 
Furthermore, Knill shows that the extensions that we treat in this paper to non-nearest 
neighbor interactions and conditionally applied gates are in fact included in Valiant's class 
of gates, albeit in a nonconstructive manner. 

In Section we establish the mapping from Valiant's gate set to a system of fermions. 
In Section |ITl| we show how the classical simulation comes about when we restrict ourselves 
to quadratic interactions that preserve the fermion number and in Section |^ we handle 
the general case of noninteracting fermions. Finally, in Section [V| we show how classically 
conditioned gate operations can likewise by simulated with our methods. 



II. NONINTERACTING FERMIONS 



Let us first state the main theorem of Ref. we will give a slightly restricted version of 
the theorem that does not include the extra freedom of gate choice on the first two qubits 
nor the possibility to do non-unitary gates: 

Theorem 1 (Valiant [^) Let M be the unitary transformation representing a quantum 
circuit on n qubits that consists of 2-qubit gates U on qubits Xi and Xi+i, i = 0, . . . ,n — 1, 
where e^'^U is of the form 



^ Ul^ ul^^ 
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(1) 



where and U"^ are arbitrary elements of SU{2) and is an arbitrary phase. There exist 
polynomial-time classical algorithms that evaluate (1) \{y\M\x)\^ for arbitrary bitstrings x 
and y, (2) Tr(?/*|M|a;)(a;|M^|?/*) where y* corresponds to an assignment of an arbitrary 
k-bit subset for any k, and (3) sample, given an arbitrary input string \x) , the probability 
distribution over outcomes y* of a measurement (in the computational basis) on an arbitrary 
k-bit subset of the qubits. 

Note that (3), which corresponds to the final simulation of a quantum computation, 
follows from (2) in a fairly straightforward manner (see Ref. |jl[]) whereas (2) could be strictly 
stronger than (1). 

A first observation about the class of allowed gates is that they preserve the parity of an 
input bitstring \x), which is expressed by the fact that the {|00), |11)} sector is decoupled 
from the {|01), |10)} sector. Note that the overall phase factor e**^ is irrelevant in the 
computation. 

To make contact with physical models, we write a gate U acting on nearest-neighbor 
qubits i and i + 1 as e*^. This Hamiltonian H can be written as a sum of three types of 
interactions: 

Hi = aiZi® + (3ili ® Zi+i, (2) 
H2 = a2Xi®Xi+i + (32Y,®Yi+i, (3) 

and 

i/s = a^Xi ® + [3^Y,(d X,+i . (4) 
where ttj, j3j are real, and X, y, Z are the three Pauli matrices: 
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At this point we note that the gate set in Theorem |I] seems extremely close to a universal 
set of gates. It has been proved |^ that universal quantum computation can be achieved by 
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employing only the Xy-interaction, i.e. HocX^X + Y^Y, if these gates can be applied 
on any pair of qubits. Since this form of interaction is certainly allowed in Theorem we 
conclude that the nearest-neighbor constraint is crucial in the construction. Another obser- 
vation is that adding arbitrary 1-qubit gates to this gate-set would also result in universality; 
it has been proved that universal quantum computation can be obtained with a circuit with 
arbitrary 1-qubit gates and only nearest-neighbor (1 dimensional) Xy- interactions 

Let us now consider the mapping onto a system of fermions. We can identify the n- 
bit computational basis states \x) with a state of n fermionic modes, each of which can be 
occupied, corresponding to 1, or unoccupied, corresponding to 0. We have a set of operators, 
creation operators aj and (hermitian conjugate) annihilation operators a, associated with 
each mode i which obey the anticommutation rules 

{ttj, Uj} = UiUj + UjUi = 0, {a|, aj} = 0, {oj, a]} = SijI. (6) 

The annihilation and creation operators act on computational basis states in the following 
manner, consistent with the anticommutation relations: 

(^i l-^Oi . . . Xj, . . . , Xii—l'} ^Xi,l(^ m—O |xo, • • • , 3^1, • • • , 2-ri— l) , (7) 

and 

aJ|xo, ...Xi,.. .,Xn-l) = 5x,,oe*''®'"=o^'|xo, ...,Xi,.. .,Xn-l). (8) 

Given these definitions, we can transform the Pauli operators in Eqs. (0-^ to creation and 
annihilation operators of fermions by a Jordan- Wigner transformation 0,^. This is done 
by first defining the operators af = ^{Xi ± zYi) that relate to the annihilation and creation 
operators as 

a+ = e^" S-i a] , a J = e'^ . (9) 

With these rules, the three types of interactions Hi,H2 and H3 can be rewritten as (we 
omit terms that are proportional to I since they will only add irrelevant phase factors to 
the quantum state of the computer): 
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and 



Hi = 2aialai + 2/5ia|+iai+i, (10) 
H2 = 02(4 - + aj+i) - P2ial + ai)(aj+i - Oj+i), (11) 



Thus we see that the total Hamiltonian H = Hi + H2 + is a sum of nearest-neighbor 
fermionic interactions that are quadratic in the fermion creation and annihilation operators, 
i.e., we can obtain any real linear combination of the Hermitian operators ajaj, aj_,_^aj+i, 
a|aj+i — ajfll+i, i{alai^i + aial_^_i), alal_^_i — ajOj+i, and i(a|a|^_]^ + ajaj+i). Fermionic systems 
that evolve according to such a quadratic Hamiltonian are referred to as 'noninteracting', 
because a canonical transformation (change of basis) exists that brings the Hamiltonian into 
a standard form involving a sum of terms each of which acts only on a single mode. 

We note that if the initial gate set in terms of Pauli matrices did not have the nearest- 
neighbor restriction, then the corresponding fermion interaction would not have been 
quadratic: this is due to the nonlocal 'sign' part in the Jordan- Wigner transformation, 
Eq. (P). It has been found that these nonlocal signs are not a problem when one consid- 
ers the question of simulating the dynamics of fermionic systems on a quantum computer: 
fermions dynamics can be simulated efficiently on a quantum computer, see Refs. Q, 



and |Tl|. Furthermore, it has been shown in Ref. |Ty] that universal quantum computation 
can be obtained by fermionic interactions that include Hamiltonians that are quartic in the 
annihilation and creation operators. Terms with an odd number of fermion operators are 
unphysical (they could transform an isolated fermion into an isolated boson), but they have 
some interesting mathematical features, see the discussion. Section 

III. PRESERVING THE NUMBER OF FERMIONS 



Before we discuss how a fermionic circuit involving Hi, H2 and H-^ can be simulated 
classically, we show how this simulation is done in the more restricted case when the gates 
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preserve the number of fermions. Thus we consider a circuit on n fermionic modes where 
each elementary gate U corresponds to an interaction between modes i and j, and can be 
written as f/ = exp{iHg) where the gate Hamiltonian is written generally as 

Hg = hiia\ai + bjja'jaj + bijajaj + b*ja^jai. (13) 

Note that the coefficients baf3 form a 2 x 2 Hermitian matrix; we will consider these coefficients 
to be part of an n x n matrix b, which is only non-zero for matrix elements involving modes 
i and j. Here and later in Section ^ we impose no restriction that i and j be nearest- 
neighbor modes, unlike the case that Valiant introduced. We will abbreviate the vacuum 
state |00 . . . 0) as |0). Let U = f/poiy(n) • • • U2U1 be a sequence of 2-qubit gates representing 
the quantum circuit. We consider 

f/a||0) = UalU^U\0) = UalU^O), (14) 

since U\0) = |0) due to fermion number preservation. U acts by conjugation as 

UalU^ = J2V,„,al. (15) 

m 

When U corresponds to a gate operation as in Eq. ([T3|) , the matrix V is given by \^ = 
exp(ib). This result is proved by making a canonical transformation that diagonalizes b. 
By group composition, the matrix V for the entire circuit U is given by matrix multiplication 
of the Vs for each gate. This evaluation of V is polynomially efficient in n if the circuit 
contains poly(r;,) gates (In fact, we could replace the individual 2-qubit gates by an arbitrary 
quadratic fermion-number preserving Hamiltonian and the matrix V of the total circuit could 
still be evaluated efficiently). 

We will ffist show how to evaluate efficiently the matrix element {y\U\x) where \x) and 
\y) are arbitrary input and output bitstrings. Since U preserves the number of fermions, 
{y\U\x) = if X and y have different Hamming weight. Let U act on a state |x) with k 
fermions in certain positions: 

U\x)=Ualal...al\0), (16) 
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where zi < 22 < • • • < "^fc by convention. Using Eqs. ([T^JT5|) , we can write 

u\x) = Vh,PlV^2,P2 ■ ■ ■ Vik.pAA^ ■ ■ ■ (1^) 

The output state equals {y\ = {0\ai^ ■ ■ - dh, where h < I2 < ■ ■ ■ < h- Using the anticom- 
mutation rules Eq. we see that contributions to the inner product {y\U\x) only arise 
when pi . . .pk is some permutation vr of the indices li . . . Ik- Furthermore, we get an overall 
sign for every such term corresponding to the number of interchanges of creation operators 
that we have to perform in order to rewrite the state o-l-ih) ■ ■ ■ '^!r(«fe)|0) '^h • • • Af- 
ter this reordering no more sign changes will take place, since (0|aij, . . . ai^al^ . . . a] JO) = 
{0\aii,a]^ . . .a/jaJjO) = 1. Thus 

{y\U\x) = sign{n)Vi,,^(^i,)Vi^,^(^i^) . . . (18) 

If V is defined as the matrix V where we have selected rows ii, . . . ,ik and columns /i, . . . , h, 
then we see that {y\U\x) = det(V^). The determinant of a kxk matrix, k < n, is computable 
in polynomial time in n. 

A. Simulating Measurements 

Next we consider how to simulate classically the outcomes of measurements on arbitrary 
subsets of qubits at the end of the computation. We will show how to calculate the prob- 
ability that a certain subset of qubits is in a particular state y* (item (2) in Theorem p. 
With those probabilities in hand, one can sample the probability distribution as given by 
quantum mechanics (item (3) of Theorem |ID. 

The Hermitian operator afoj counts the number of fermions in mode i. Thus its expec- 
tation value with respect to a density matrix p, Tr ajaip, is the probability that mode i is in 
state Similarly, the expectation value of a^al = I — a\ai is the probability that mode i is 
in state |0). So in order to evaluate the probability that, given an input state a certain 
subset of k modes is in state \y*) we calculate 
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p{y*\x) = Tiaj^a]^ . . . a]^aj^U\x){x\U^ = {x\ U^aj.a]^ . . . a]^aj^U \x), (19) 

with ji 7^ j2 7^ • • • 7^ jk and we use a|aj or aja| when y* = 1 or respectively. (Eq. ([I9D 
illustrates a case where y*_^ = and y*^ = 1.) Again we write \x) = aj^ . . . ctpjO). We have 
to evaluate the expression 

P{y \^) = ^ ^ji,mi^ni,ji ■ ■ ■ ^jk,mk^nk,jk 

mi,ni,...,m,i^,nk 

{0\ap^ . ..ap,{an,al^ ■ . .a|„^a„jaj^ . . .aJjO). (20) 

We will invoke Wick's theorem [|12| (described in quantum many-body or quantum field 
theory textbooks such as Ref. [|13[) to rewrite this formula^ Wick's theorem states that we 
can rewrite a string of annihilation and creation operators Ai ... An as 

\n/2\ 

A^...A^ = : A^...An: + J2 Ck, (21) 

fc=i 

with 

Ci = : ...An: + : AMA3 . . . v4„ : + . . . 

C2 = : AA2 A^A4 A5...An: + : A^A^A^ A5 . . . A^ : + . . . 

etc. (22) 

Here : Ai . . . An : denotes the so called normal ordered form of the sequence of operators 
Ai . . . An. : Ai ... An : is equal to the reordered sequence of operators . . . where 
all the creation operators are moved to the left (but not reordered amongst each other), 
and the quantity is negated when the number of interchanges of creation and annihilation 



We actually use Wick's theorem for ordinary operator products rather than for time-ordered 
operator products, which is the case analyzed in most standard treatments. The version of the 
theorem we are using is excellently set out in T. D. Crawford and H. F. Schaefer III, "An In- 
troduction to Coupled Cluster Theory for Computational Chemists", Reviews in Computational 
Chemistry 14, 33-136 (2000); see also |http://zopyros. ccqc.uga.edu/lec_top/cc/html/nodel3.html . 



9 



operators to achieve this form is odd. The object AiAj is called a contraction and is defined 

as 

I — I 

The terms Ck in Eq. (|2ll) are each a sum over every possible choice of k contractions in the 
normal ordered product. 

From the anticommutation rules for creation and annihilation operators, it follows that 

rn rn t i i „ i 5. 

alttj = alttj = ttittj = U, ttittj = Oij. (z4j 

The normal-ordered form is extremely convenient when evaluating an object such as 
(0|Ai . . . y4„|0) since the vacuum expectation value of any normal-ordered sequence of opera- 
tors vanishes (0| : Ai-^ . . . Ai^ : |0) = 0. Therefore, when we evaluate the vacuum expectation 
value of Eq. ([2^), the only terms that remain come from C\n/2\ even), in which every 
operator Ai is contracted (or matched) with another operator Aj. The last step is to bring 
the fully contracted terms to a form in which contracted operators are adjacent, that is, we 
have 



(0| : A^A^A^A^A^ ...A^... A^^^A^ : |0) 



I 1 



sign(7r)(0| : A^A^^ A^A^ A^A^ . . . A^-xA^ : |0) = sign(7r) A^A^ A^A^ A^A^ . . . A^.^A^, (25) 

where sign(7r) is —1 (1) when the number of crossings of the contraction lines is odd (even). 
Evidently, what emerges is the Pfaffian Pf(M) of M{i,j), an n x n antisymmetric matrix 
(i.e. M{i,j) = —M{j, i)). The Pfaffian Pf(M) is when n if odd, and for even n it is defined 
asQ 



^For an discussion of Pfaffians via an exposition of Pfaff 's 1815 paper introducing the construction, 
see T. Muir, The Theory of Determinants in the Historical Order of Development (Vol. 1, Dover, 
1950), pp. 396-401. 
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Pf (M) = ^ sign(7r)M,(i),,(2) . . . 

(n— l),7r(n) ; 

(26) 

TT 

where the sum over vr is restricted to permutations on the indices 1, 2, ... n such that iT{2k — 
1) < 7T{2k) and 7r(l) < 7r(3) < 7r(5) . . .. Eqs. (|0|,|5],|6D tell us that 

p{y*\x)=Pf{M), (27) 

where M can be constructed from Eq. (^) and the contraction identities, Eq. (^), in the 
following manner. The matrix elements M{i,j) for 1 < i < j < 2{k + I) are obtained from 
Table 1: The indices i,j = 1, . . . ,2{k + I) are assigned to the ordered sequence of creation 
and annihilation operators in Eq. (^). To determine M{i,j) {i < j) we find what type 
of operator the indices i and j correspond to and then read off the matrix element M{i,j) 
from the table. We use unitarity of V and the contraction rules to determine each entry of 
the table. The Xs in the table indicate that these entries do not occur. 

The Pfaflian of an x antisymmetric matrix M can be computed in poly(n) time, 
since Pf(M)^ = det M. The simulation procedure that was formulated in Ref. |]1| very 
similarly relies on the evaluation of a Pfaffian. At the moment it is not clear to us how 
the representation of the quantum circuit in Ref. Jll using matchgates corresponds to the 
fermionic representation developed here. 

IV. GENERAL NONINTERACTING FERMIONS 

We are now ready to consider the classical simulation of a quantum circuit consisting of 
gates that are built from general quadratic fermionic interactions. These interactions only 
preserve the parity of the photon number. In order to deal with these general interactions, 
we transform the set of fermion annihilation and creation operators to a new set of Hermitian 
operators (associated with so called Majorana fermions |TD|,|T3): 



C2i = ai + al, C2i+i = -i{ai - al) , (28) 

where i = 0, . . . ,n — 1. The anticommutation relation for this new set of operators is 
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{ck,ci} = 2SkiL (29) 

Note that operators C2i and C2i+i are in some sense the fermionic version of conjugate vari- 
ables p and q that are obtained from hnearly combining bosonic annihilation and creation 
operators. It is clear that the Hamiltonians Hi,H2 and of Eqs. (0-^ will be quadratic 
in these new operators. Let f/ be a sequence of 2-qubit gates each composed of interactions 
that are quadratic in the operators q, i.e. each of the gates corresponds to a Hamiltonian 
H 

H = -^akiCkCi. (30) 

We again have omitted any term proportional to /. Hermiticity of H requires only that 
lm{aki) = lm.{aik)- It is conventional to choose the a matrix to be real and antisymmetric, 
so that i(x is a Hermitian n x n matrix. For the interactions we have introduced, the a 
matrix will only be nonzero in a 4 x 4 subblock, but this restriction is not necessary for the 
following procedure to work. Similar to the number conserving case, a sequence of gates 
U = t/poiy(n) . . .U2U1 acts by conjugation as 

Uc,U^ = Y.R,jC,, (31) 
j 

where R G S0{2n). We will establish this important result by explicitly computing the 
matrix R for a single gate U = e^^ . The result is not so well known as for the number 
conserving case (although it has been mentioned in [|1^), so we will give some of the details 



of the derivation. We follow the notation of [|1^. First, the Hamiltonian of Eq. (pOD, with 



CK chosen to be a real antisymmetric matrix, can be brought into canonical form 



h=-,y: ^m- (32) 



2. to 



h' and h" are given by the real orthogonal transformation^ 



^See R. A. Horn and C. R. Johnson, Matrix Analysis (Cambridge, 1985), p. 82, Theorem 2.3.4. 
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h" 



h" 



w 



Co 
Cl 

C2n-2 
y C2n-1 y 



(33) 



The 2n x 2?7, orthogonal matrix W diagonahzes cx into 2x2 blocks: 



WaW^ 



eo 
-eo 



\ 



V 



(34) 



e„-i 

-Cn-l 

The 6's and 6"s have the same anticommutation relations as the original Majorana fermion 
operators. Note that ie^ are the eigenvalues of the matrix ia.. We now write Eq. ( pi]) using 
the canonical transformation: 

^Qf/t = Y^eM-\T.^'n.h'j>i){w,,,h'^ + w,,^,,h';)eM\T.^n.h'j:^). (35) 

j m m 

Because the Hamiltonian in this canonical form is a sum of commuting terms, the exponen- 
tials here can be factorized; the factors with m ^ j commute through and disappear, and 
we obtain 



f/Qf/t = Y^eM-\^M){W,,,h'^ + W^2,+i.6;')exp(le,-6;.6;0- 
The remaining exponential factors can be expanded and simplified: 



(36) 



k — k — 



00 / /9^2fc 00 / /9\2fc+l 



- (2A:+1)! 



cos(ej/2) + 6'6;'sin(e,/2). 



(37) 



Plugging this form into Eq. ( PB] ) and simplifying, we obtain Eq. (PT|), where R is given by 
R = W'^MW, with the matrix M having the 2x2 block form 
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M 



cos eo — sm eo 
sin eo cos eo 



V 



(38) 



cose„_i -smen_i 
sin en-i cos e„_i 

As before, given a quantum circuit with poly(n) two-mode noninteracting fermion gates, 
that is, involving four Majorana fermions per gate, we can construct the total 2n x 2n 
matrix R in polynomial time by straightforward matrix multiplication of the individual R 
matrices corresponding to the gates. 

We again consider the probabilities with which certain measurement outcomes are ob- 
tained, i.e. p{ii*\x) in Eq. (|T^, and show that as before these quantities are equal to the 
Pfaffian of some antisymmetric matrix. 

As before, we consider an input state |a;) = o)p^ ■ ■ .aj,jO), which we will now write as 
C2pi • • • C2pjO) with pi < p2 < . . . pi- Thus we would like to evaluate 



p{y*\x) = (0| C2pi ■ ■ ■ C2pi 



(39) 



The pattern of a and is again determined by the state |?/*). We need a formula for how 
our general (non-number conserving) U acts by conjugation on the creation and destruction 
operators. We can use Eq. (|3T|): 



U^a,U = ^f/t(c2. + ^C2^^l)U = ^ Y^i^h + ^^^.+1,.)^. = E T^,c,, (40) 



and similarly 



U^alU = Y.T*,c,. 



(41) 



This defines the n x 2n matrix T. We then obtain for the measurement probability: 



p{y*\x) 



E 



Ji.mi-' jl,ni 



mi,ni,...,mfc,nfe 



(0|c2pj . . . C2pj CjT^jCjij . . . Cnf.CYnf.C2pi ■ ■ ■ C2p; |0). 



(42) 
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Again we can use Wick's theorem to evaluate the vacuum matrix element. This is done by 
writing the Majorana operators in terms of the fermion creation and annihilation operators. 
Expanding gives a large number of terms, to each of which Wick's theorem applies. Each 
term normal orders differently, but in every case only the fully contracted terms survive. 
All of these fully contracted terms are generated by contractions directly over the Majorana 
operators, defined by linear extension: 



C2iC2j+l 



-t[aiaj 



I ^2 



a j cij I 



i5. 



(43) 



And similarly 



I I 

C2i+lC2j 



—i6. 



n _ I I 

C2iC2j — C2i+lC2j+l 



(44) 



Here we have used Eq. (p^. Then, the vacuum expectation value is written as the sum of 
all fully contracted expressions over the Majorana operators, with the usual fermionic sign. 
Thus, we can say that Wick's theorem applies in the same way to the Majorana fermion 
operators as it does to ordinary fermion creation and annihilation operators; we emphasize 
that this is only true for the vacuum expectation value, it is not true as an operator identity 
(normal ordering is not defined for the Majorana operators). 

We can summarize these contraction rules by writing CjCj = Hij where H is a. 2n x 2n 
Hermitian matrix consisting of 2 x 2 blocks: 



H 



1 i 
-i 1 



v 



1 i 
-i 1 



(45) 



Applying Wick's theorem again leads to a Pfaffian expression. In Table II we give entries 
that permit the 2{l + k) x 2{l + k) matrix to be constructed such that p{y*\x) = Pf(A^). 
Again, the entire evaluation is clearly doable in polynomial time. 
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V. INTERMEDIATE MEASUREMENTS AND CLASSICALLY CONDITIONED 

OPERATIONS 

We now extend our quantum circuit of noninteracting fermions by allowing intermediate 
complete von Neumann measurements in the computational basis on subset of qubits, which 
then determine the subsequent choices of unitary gates and measurements on the remaining 
qubits. We will show here that a fermionic circuit with these resources can still be simulated 
efficiently with a classical algorithm. Care has to be taken in specifying which intermediate 
and final measurements are allowed in our model; we restrict ourselves to complete von 
Neumann measurement in the computational basis, i.e. the outcomes of the measurement 
are either 'no fermion present in this mode' or 'one fermion present in this mode'. A lot 
of added power can be hidden in the kind of measurements that one is allowed to do; for 
example it has been shown in Ref. [0 that universal quantum computation can be achieved 



by noninteracting fermion gates plus a nondestructive eigenvalue measurement of the quartic 
operator cjCkCrCs- 

A general quantum circuit employing our set of resources is depicted in Fig. 0. Every 
time a measurement is made on a subset of qubits, these qubits are no longer used in 
any later steps of the computation. Our classical simulation will be constructed in the 
following manner. Measurements Mi, . . . , M^, on subsets Si, . . . , will take place at 'times' 
The total unitary evolution until the measurement Mi is denoted as Ui, the 
conditional unitary evolution between M^ and M^+i is denoted as Uk+i{yl, yl, . . . , yl) where 
the labels yl,y2, ■ ■ ■ ,yk correspond to the outcomes of the measurements Mi, . . . , M^. The 
choice of measurements themselves may depend on earlier measurement outcomes, i.e. Mi = 
Mi{yl, . . . , yi_i). Even though the later time-evolution operators will not act on the qubits 
that are already measured, we keep the dimension of these matrices the same as the initial 
matrix Ui, i.e. these are 2" x 2" matrices when the total number of fermionic modes is n. 

We can calculate the probability that at time ti subset 5*1 is in the state \yl) (and sample 
from this probability distribution) by the methods that we have developed in Sections fT| 
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and If the quantum measurement gives \yl), then the remaining qubits are in the state 



^ UM)PrUi\x){x\ulPrU2iylV 

TTPy,U,\x){x\Ul ' ^ ^ 

where the projector Py* is of the form aj^a]^ . . . aj^^ I'^iisii where G 5*1, and 

whether the factor djiOj- or Oj.aj. appears depends on whether {yl)j^ is or 1 respectively. 
Let us assume that we have sampled the measurement probability distribution at time ti 
and have found a particular outcome yl. To sample from the probability distribution of 
measurement M2, we will have to be able to evaluate 

p{y;\ylx) = TTPy*^P2. (47) 

Py*, like Py*, is again a product of creation and annihilation operators, e.g., 
al^Oj^ . . . a||^^|aj|g^| where ^i, . . . ,^1521 ^ 'S'2 and the pattern of creation and annihilation op- 



erators depends on the bits of 1/2 • The denominator in Eq. p6| ) is already determined when 
simulating the first measurement, so we will focus on calculating 

Tr Py,UM)PylUi\x){x\UlPy*U2{yiy = {xlUlPy^UMV Py*UM)PylUi\x). (48) 



This equation has basically the same form as Eq. (PDf), except that (1) we have more anni- 
hilation and creation operators and (2) we conjugate different sets of operators by different 
unitary matrices. The important fact here is that we can again express the probability as 
the Pfaffian of some antisymmetric matrix. Let us see how we construct this matrix. 

At this point we simplify the notation somewhat. Let U2{yl) = U2, Py* = Pi and 



Py* = P2. We put in UiUl and U2U2 terms in the appropriate places in Eq. (^Sf), so that 



operators in the first (from the left) Pi get conjugated by Ul, operators in P2 get conjugated 
by UI2 = U1U2 and the last Pi gets conjugated by Ul again. Let T^, k = 1, 12, be defined 
by UlaiUk = J^jTjjCj. We obtain: 

Piyh y2\^) = ^ii,ai^/i*fei • • • ^j|Sj|,fe|s^|^j|Si|,a|5^| 

ai,i)l,/l,fll,....a|g^|,b|5;^|,/|g^|,S|3^ 

di,ei,...,d|52|,e|52 
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rpXI * r-p\1 rpVl * rpl2 rjil rjil * rpl * rjil f AQ\ 

-^iuei-^iudi ■ ■ ■-^i\S2\^e^S2\ i\S2\As2\ hJl h,9l ' ' ' ^ 3\S^\,9\S^\ 3\S^\J\S^\ > 
(0|c2p, . . .C2pi (CaiCfoi . . . Cfe|g_^|Ca|g_^|)(CeiCdi . . . Ce^g^^Cd^^^^){Cf^Cg^ . . ■ Cg^s^^Cf^g^^)c2p^ ■ ■ .C2pj0). 

In Table III we show how to construct the matrix O of dimension 2(/ + 2|S'i| + 15*21) for 
which p{yliy2\x) = Pf(0). The notation C(^a/f)a indicates that the c-operator can be either 
a Ca^ or a c/^; the reason is that these operators have identical T prefactors. 

It is clear that we can extend this procedure to the case of a circuit that contains 
k = poly(r;,) instances of measurements on subsets that determine the next choice of unitary 
evolution. In general when we express a probability such as Eq. (4S), we see that Pi gets 
conjugated by Ui, P2 by Uu, P3 by f/123 and Pk by Ui2...k, the total unitary evolution. 
When we write p{yl, 1/2 ■ ■ - Vkl^) = Pf(^)) the dimension of the matrix X is 2(/ + 15^1 + 
2I^i=i I'S'jl). The entries of this matrix can be determined by calculating particular matrix 
elements (specified by the measured sets of qubits) of at most {2k + 1)^ matrices of the form 
T^HT^'^ etc. where i and j are labels which can be 1, 12, 123, . . . , 123 . . .k. 

Let us summarize the simulation algorithm: 

Classical simulation of a quantum circuit with noninteracting fermions and 
fermion counting measurements, see Fig. [l|: 



1. Compute the n x 2n matrix corresponding to Ul, Eq. (^Ol). 

2. Simulate measurement Mi. sample from the probability distribution p{yl\x) using the 
measurement theorem in Ref. and the fact that p{y*\x) = Pf(Oi) = Vdet Oi where 
Oi is a 2(Z + A;) x 2{l + k) matrix with k equal to the Hamming weight of input string 
X and / equal to the number of bits in y^. 

3. Let yl be the outcome of this measurement Mi, and let U2 be the corresponding 
unitary evolution. Compute T^^ corresponding to UiU2- 
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4. Simulate measurement M2: sample from the probability distribution p{y2\yi,x) — 
^^p{y*fl)'' where we use the fact that we can evaluate p{yl, 1/2 1 a;) = Pf (02)- O2 depends 
on T^^ and as in Table 111. 

5. Let 7/2 be the outcome of the measurement M2 and let U3 be the corresponding unitary 
evolution, possibly also depending on the first outcome y^. Calculate T^^^ correspond- 
ing to uMul 

6. Simulate measurement M3: sample from the probability distribution p{ys\yi,y2,x) — 

^^p{y*'y*\^)'' ' "^^^^^ '^^ ^^^^ ^^^^ evaluate p{yl,y2,y3\x) — Pf(03). O3 

depends on T^^^, T^^ and T\ 

7. Repeat steps 5 and 6 for the subsequent evolutions U4, . . . , Uk, finding expressions for 
21123... fc^ and finally simulate the last measurement by sampling from the distribu- 
^ionp{yM,yl...,yU,x) = 0g^y 

It is evident that this procedure is polynomial when the number of stages k of the 
compute/measure procedure of Fig. 1 is poly(n): the largest matrix whose Pfaffian must be 
computed has dimension bounded by Akn. 



VI. DISCUSSION 

The present work opens a set of very interesting questions concerning the boundary 
between classical and quantum computation. For fermionic quantum circuits we may ask, 
what is the effect of adding circuit elements beyond those considered above (those associated 
with a noninteracting fermion model)? Three outcomes are possible: 1) the circuit can 
perform universal quantum computation, 2) the circuit remains efficiently simulatable by 
a classical computation, or 3) some intermediate case. For example, one could explore the 
effect of adding (unphysical) hnear terms to the gate Hamiltonians. These terms a, + aj 
and i{ai — a|) will be somewhat similar, but not identical to 1-qubit gates in terms of Pauli- 
matrices; they are nonlocal gates as can be seen from the Jordan- Wigner transformation, Eq. 
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(^). It can in fact be shown that these hnear interactions can be incorporated in purely 
quadratic fermion interactions by adding a new fermionic mode, which we may label " — 1", 
and changing the linear interactions on, say, mode i to quadratic interactions between mode 
i and mode —1. Another line of investigation could be into some known physical models 
(e.g., the Anderson model and the Kondo model |TB[) that involve more general fermionic 
interactions at a single site (meaning 2 fermionic modes), many of whose properties are 
computable. 

Valiant's work shows that some terms added to the fermion model (or some gates added 
to the circuit model) still result in a classically simulatable system; these new terms are to 
act on the first two fermionic modes only. These new gates (if they exist) do not preserve 
the parity of the number of fermions, and thus involve either linear or cubic terms in the 
annihilation and creation operators. (This follows from the fact that any gate that obeys 
the 5 matchgate identities in Ref. |jl|] and preserves the parity of the number of fermions, 
is automatically quadratic in the number of fermions.) See Ref. |^ for a more extensive 
treatment of these additional gates. 

The case of adding power by intermediate measurements is also quite interesting: be- 
tween the case of complete von Neumann measurements that are classically simulatable, 
and the quartic-operator basis measurements of Bravyi and Kitaev ||10|, which give uni- 
versal quantum computation, there are many possible POVM measurement scenarios that 
have not been analyzed. We are hopeful that further analysis will be able to identify more 
scenarios as definitely classical or definitely quantum. 

This classical-quantum boundary is remarkably different for fermionic and bosonic sys- 
tems. A model quadratic in bosonic operators describes the "linear optics" scenario of quan- 
tum computation; formally, this model only differs from the noninteracting-fermion model 
in that bosons are characterized by commutation rather than anticommutation relations: 

[oj, Qj] = ttittj — ajai = 0, [a|, a]] = 0, [oj, a]] = 6ijl. (50) 

Nevertheless, the exact parallel to the model we analyzed above, in which quadratic Hamil- 
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tonians can be interspersed with complete von Neumann measurements, is fully "quantum" 
in the bosonic case despite being "classical" for fermions. 

We would like to emphasize again that this quantum/classical distinction may not be 
perfectly sharp; being able to efficiently compute some properties of a circuit classically 
does not mean that every aspect of the quantum dynamics of this circuit are also efficiently 



computable. This is shown when we try to carry out the analysis in Section |T| and Section 
for bosons, to see where the parallels between the two cases break down. In fact, one 
can go surprisingly far before any differences appear. An equation of the form of Eq. (p!5D 
still applies for bosons, where V can be be shown to be any element in U{n) [|l^. That 
is, a canonical transformation exists in both the boson and fermion case which permits a 
representation of the system as a set of noninteracting particles. 

It might have been thought that this easy diagonalization leads to an easy computation 
of all possible properties, but this does not appear to be the case. An essential difference 
is that no sign changes occur when we interchange the bosonic creation operators amongst 
each other. This causes an expression such as {y\U\x) in Section ^ to be equivalent to the 
permanent of some matrix, if we analyze the case when |x) and \y) both contain no more 
than one boson per mode. The permanent is a much harder object to calculate exactly than 
the determinant of the fermion case, in fact this has been proved to be an ^^P-complete 
problem ||T^. Our methods will therefore fail to evaluate these bosonic matrix elements 
efficiently. 

In summary, our results on the classical simulatability of noninteracting fermions leave 
some interesting questions unanswered about the computational power of various physical 
models of quantum computation. While Valiant's analysis turns out to conform largely to a 
"known" area of physics, his work shows that mathematical approaches to these problems 
are possible that have never been envisioned in many-body physics. 
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FIG. 1. A quantum circuit with classically conditioned gates. 
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TABLE I. The matrix M{i,j) for i < j. 





3 








C2p0 
















{T*H)j^,2p0 


C2p„ 




{HT'^)2pa,j0 


Sa,f3 



TABLE IL The matrix N{i,j) for i < j. 
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TABLE in. The matrix 0{i,j) for i < j. 
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